Next‐generation matrices for marine metapopulations: The case of sea lice on salmon farms

Abstract Classifying habitat patches as sources or sinks and determining metapopulation persistence requires coupling connectivity between habitat patches with local demographic rates. While methods to calculate sources, sinks, and metapopulation persistence exist for discrete‐time models, there is no method that is consistent across modeling frameworks. In this paper, we show how next‐generation matrices, originally popularized in epidemiology to calculate new infections after one generation, can be used in an ecological context to calculate sources and sinks as well as metapopulation persistence in marine metapopulations. To demonstrate the utility of the method, we construct a next‐generation matrix for a network of sea lice populations on salmon farms in the Broughton Archipelago, BC, an intensive salmon farming region on the west coast of Canada where certain salmon farms are currently being removed under an agreement between local First Nations and the provincial government. The column sums of the next‐generation matrix can determine if a habitat patch is a source or a sink and the spectral radius of the next‐generation matrix can determine the persistence of the metapopulation. With respect to salmon farms in the Broughton Archipelago, we identify the salmon farms which are acting as the largest sources of sea lice and show that in this region the most productive sea lice populations are also the most connected. The farms which are the largest sources of sea lice have not yet been removed from the Broughton Archipelago, and warming temperatures could lead to increased sea louse growth. Calculating sources, sinks, and persistence in marine metapopulations using the next‐generation matrix is biologically intuitive, mathematically equivalent to previous methods, and consistent across different modeling frameworks.

popularized in epidemiology to calculate new infections after one generation, can be used in an ecological context to calculate sources and sinks as well as metapopulation persistence in marine metapopulations. To demonstrate the utility of the method, we construct a next-generation matrix for a network of sea lice populations on salmon farms in the Broughton Archipelago, BC, an intensive salmon farming region on the west coast of Canada where certain salmon farms are currently being removed under an agreement between local First Nations and the provincial government. The column sums of the next-generation matrix can determine if a habitat patch is a source or a sink and the spectral radius of the next-generation matrix can determine the persistence of the metapopulation. With respect to salmon farms in the Broughton Archipelago, we identify the salmon farms which are acting as the largest sources of sea lice and show that in this region the most productive sea lice populations are also the most connected. The farms which are the largest sources of sea lice have not yet been removed from the Broughton Archipelago, and warming temperatures could lead to increased sea louse growth. Calculating sources, sinks, and persistence in marine metapopulations using the next-generation matrix is biologically intuitive, mathematically equivalent to previous methods, and consistent across different modeling frameworks.

K E Y W O R D S
marine systems, metapopulation, next-generation matrix, salmon farms, sea lice, source-sink dynamics

T A X O N O M Y C L A S S I F I C A T I O N
Conservation Ecology, Population Ecology, Theorectical Ecology

| INTRODUC TI ON
Metapopulations consist of subpopulations located on isolated habitat patches that are connected via dispersal (Hanski, 1998;Kritzer & Sale, 2004;Levins, 1969). In most benthic marine species, this dispersal comes from the pelagic larval stage (Cowen & Sponaugle, 2009). Larvae disperse between and then settle on habitat patches, and once settled the remaining stages are sedentary and remain confined to a specific habitat patch. In marine systems, the metapopulation concept, where subpopulations are connected but have their own demographic rates, has been used in the spatial planning of marine protected areas (MPA) and the siting of marine reserves (White et al., 2013).
In a metapopulation framework, habitat patches are often classified into sources and sinks based on how the subpopulations on these patches contribute to the overall metapopulation. The source-sink classification of habitat patches was first described concretely in a terrestrial context by Pulliam (1988), where habitat patches were classified as sources if the local subpopulations could persist in isolation and sink if they could not. However, this classification ignores the effect of dispersal, which is especially critical in marine metapopulations, and so Runge et al. (2006) and Figueira and Crowder (2006) updated the classification of source and sink patches to include both the local productivity of a patch, as well as the ability to disperse away from the habitat patch.
Under this new classification, a source patch is a patch on which an adult will more than self-replace over the entire metapopulation, and a sink is a patch on which an adult will not. Self-replacement need not occur on the same habitat patch as the adult originated, and thus under this classification, a source patch may not be able to persist in isolation.
Due to the large-scale larval dispersal that occurs in many marine species, it is common in marine metapopulations for source patches to not be able to persist in isolation (White et al., 2019), thus preserving a persistent metapopulation, especially in the context of MPAs, often requires more than simply preserving source patches. To maintain persistent marine metapopulations it is necessary to preserve sufficient larval exchange between closed loops of habitat patches so that an average adult can eventually self-replace over multiple generations (Burgess et al., 2014;Hastings & Botsford, 2006). If sources patches self-recruit enough larvae to persist in isolation then this can be accomplished by preserving only source patches. However, if they do not it may require preserving both source and sink patches, as the sink patches may provide sufficient larval exchange back to the source patches to create a closed loop of habitat patches over which an adult can self-replace. Evaluating the persistence of marine metapopulations is therefore difficult as it requires accurate measures of larval connectivity between habitat patches, as well as accurate local demographic rates of adult stages on each habitat patch (Burgess et al., 2014). Despite the difficulties, evaluating metapopulation persistence is critical in designing marine reserves in which the protected habitat patches can persist even when outside patches are exploited (Costello et al., 2010;Garavelli et al., 2018;Puckett & Eggleston, 2016;White et al., 2010White et al., , 2013. Next-generation matrices are a useful tool that can both be used to evaluate the persistence of a metapopulation as well as identify the contribution of local habitat patches under a variety of modeling frameworks. Originally popularized in epidemiology (Diekmann et al., 1990), they have since been used in ecology to calculate the basic reproduction number, R 0 , of structured population models (Cushing & Yicang, 1994;Li & Schneider, 2002).
Next-generation matrices convert structured population models into generational time, so that the entries of the next-generation matrix are the number of new individuals produced in each age, stage, or patch, depending on how the model is structured. The individual contribution of habitat patches or evaluation of metapopulation persistence can therefore be measured for different model structures (discrete time, continuous time, etc.) under the same framework of the next-generation matrix. The column sums can be used to measure the contribution of each habitat patch over a generation and the magnitude of the dominant eigenvalue can be used to evaluate metapopulation persistence if the nextgeneration matrix is primitive (otherwise the spectral radius can be used to determine persistence if the matrix is non-negative and irreducible; Harrington & Lewis, 2020). Next-generation operators have also been used in ecology to calculate source and sink regions in heterogeneous environments (Harrington & Lewis, 2020;Huang et al., 2016;Krkošek & Lewis, 2010;Mckenzie et al., 2012), and to evaluate the level of control required to suppress invasive species (Lewis et al., 2019).
In this paper, we focus on using the next-generation matrix to evaluate local patch contribution and metapopulation persistence in marine metapopulations, but the framework used here is also applicable to many other birth-jump metapopulations (Hillen et al., 2015), where there is a single juvenile stage which can disperse between habitat patches and the remaining stages remain on a single habitat patch. Examples of non-marine species that exhibit this structure include plant species where seeds are carried between habitat patches (Husband & Barrett, 1996), or insect species with a single large dispersal event such as the spruce budworm (Williams & Liebhold, 2000) or mountain pine beetle (Safranyik & Carroll, 2007).
In fact, the next-generation matrix approach can even be extended to metapopulations in which adults also disperse, though the calculations become more complicated and so here we focus on species with a single dispersing stage.
Specifically, to demonstrate the utility of the method, we use the next-generation matrix to calculate the contribution of a single salmon farm to the spread of sea lice in a salmon farming region on the west coast of British Columbia. Sea lice (Lepeophtherius salmonis) are a parasitic marine copepod that feed on the epidermal tissues, muscles, and blood of salmon (Costello, 2006). With a free-living larval stage, they can disperse tens of kilometers, spreading between salmon farms in a region and between wild and farmed salmon (Krkošek et al., 2006;Peacock et al., 2020;Stucchi et al., 2011).
Lesions and stress from high sea lice infestation make adult salmon more susceptible to secondary infections, leading to large economic consequences for the salmon farming industry (Abolofia et al., 2017). On wild juvenile salmon, infestation with sea lice can lead to mortality and elevated exposure to sea lice from salmon farms can contribute to population-level declines in pink salmon (Krkošek et al., 2007). In the context of sea lice on salmon farms, we are not concerned with preserving a persistent metapopulation of sea lice parasites, but instead, we use the next-generation matrix to calculate the relative contribution of each salmon farm and evaluate the effect of environmental variables on the overall growth of the sea louse metapopulation.
The specific salmon farming region that we focus on to calculate farm contribution is the Broughton Archipelago. The Broughton Archipelago is located on the west coast of Canada, between Vancouver Island and the mainland of British Columbia, and has been central in evaluating the effect of sea lice from salmon farms on wild salmon (Brooks, 2005;Krkošek et al., 2005Krkošek et al., , 2007Krkošek et al., , 2011Marty et al., 2010). The area has historically had around 20 active salmon farms (Foreman et al., 2015), though currently certain farms are being removed from this region in an agreement between the government of British Columbia and the Kwikwasut'inuxw Haxwa'mis, Namgis, and Mamalilikulla First Nations (Brownsey & Chamberlain, 2018). After 2023 many of the remaining farms must be approved by both the local First Nations and the government in order to continue to operate and thus determining the farms which are acting as the largest sources of sea lice is critical during this transition period.
The paper is structured as follows. First, we demonstrate how to use the next-generation matrix to calculate the contribution of local habitat patches to the metapopulation and evaluate metapopulation persistence. Next, we highlight how to construct the nextgeneration matrix for different types of models. Then we calculate a next-generation matrix for sea louse populations in the Broughton Archipelago to identify which salmon farms are the largest sources of sea lice in this region, evaluate the effect of the current farm removals, and investigate the effect of environmental variables on metapopulation growth. Finally, we discuss how the calculations of patch contribution and metapopulation persistence from other studies compare to the calculations using the next-generation matrix.

| MATERIAL S AND ME THODS
In this section, we present details on the construction of the nextgeneration matrix and show how it can be used to determine the contribution of the subpopulation on a single habitat patch to the metapopulation as well as determine the persistence of the metapopulation. We then present the explicit construction of the nextgeneration matrix for models with age-dependent demography and present the construction for ordinary differential equation models and discrete-time models in Appendices 1 and 2, respectively.
Finally, we detail the construction of the next-generation matrix for a system of sea lice populations on salmon farms in the Broughton Archipelago.

| Next-generation matrices
Here we construct the next generation matrix for a single species marine metapopulation with a single larval stage that can disperse between patches and where the remaining stages are confined to the habitat patch on which the larvae settle. We assume that the terminal stage is the only stage that produces new larvae. This assumption can be relaxed, though the entries of the next-generation matrix become slightly more complicated. We construct the nextgeneration matrix for models where the metapopulation is divided into l patches and the f attached stages are modeled explicitly. The larval stage is modeled implicitly, so that the birth rate into the first attached stage includes both the birth rate of larvae and the probability of larvae successfully dispersing between patches and attaching on a new patch. The lifecycle diagram for such a metapopulation is shown in Figure 1.
We use the next-generation matrix to calculate the local patch contribution to the metapopulation in the context of low population density. To calculate the next-generation matrix it is necessary to linearize a potential density-dependent model around the zero equilibrium so that the effect of density dependence at higher population sizes is ignored when calculating patch contribution. This approach of ignoring density dependence is common when determining persistence or patch contribution of marine metapopulations, as the focus is either on determining if a metapopulation can persist at all, or determining which habitat patches are acting as population sinks and which are acting as population sources (Burgess et al., 2014;Harrington & Lewis, 2020;Hastings & Botsford, 2006;Krkošek & Lewis, 2010). Alternatively, it is also useful for determining patch contribution in metapopulations of species that are being actively controlled to remain at low densities, such as sea lice on salmon farms, which is our focus in Section 2.4. Another common assumption in the theory of persistence also made here is that there is no F I G U R E 1 The lifecycle graph for two patches in a metapopulation of a species with a single larval stage that disperses between habitat patches and f sessile stages that remain on a habitat patch. The population on patch i in stage n is given by P i n . discernible Allee effect in any of the patches (Burgess et al., 2014;Hastings & Botsford, 2006;Krkošek & Lewis, 2010).
Under any model structure, the element in row i and column j of the next-generation matrix gives the total number of new individuals produced on patch i over the lifetime of one initial individual on patch j. However, how 'new' individuals are defined is subject to interpretation. Here, we construct the next-generation matrix under modeling frameworks that only consider the sessile stages explicitly, so "new" individuals will be newly attached stage 1 individuals. Then, when the number of new individuals on patch i produced from one new individual on patch j are tracked, the new individual must first survive and reproduce on patch j, before larvae disperse and arrive on patch i. In this way patch contribution, as will be calculated in Section 2.2, is primarily a function of the local patch demography which is then coupled with dispersal to other patches. Under this framework, the entries of the next-generation matrix, K, for all model structures can be given by: where the product of the first two terms is the lifetime reproduction of an individual on patch j and the third term is the probability that larvae leaving patch j successfully arrive on patch i.
However, if the next-generation matrix is constructed for models that explicitly model the larval stage, then the larval stage is often considered the first stage. The next-generation matrix will be slightly different in this case as well as the calculation of patch contribution, though the calculation of metapopulation persistence will be the same. We illustrate the differences between constructions of the next-generation matrix in the Section 4 and show how the calculation of persistence remains the same under all constructions.

| Determining patch contribution and metapopulation persistence
Here, we show how to use the next-generation matrix, K, to determine the contribution of each patch to the metapopulation and evaluate metapopulation persistence. To determine the contribution of a specific patch to the metapopulation we track the total number of new individuals produced across the metapopulation after one generation from an initial individual starting on that patch. The entries of the next-generation matrix, k ij give the number of new individuals produced on patch i over the lifetime of an initial individual on patch j. Therefore, if we define so that R j is the jth column sum of K, then R j is the total number of new individuals produced across all patches from an initial individual starting on patch j and can be used to define the contribution of patch j to the entire metapopulation.
This definition of patch contribution easily lends itself to classifying local habitat patches as population sources or sinks.
If R j < 1, then an individual on patch j cannot replace itself over the entire metapopulation, and thus patch j is defined as a sink (Harrington & Lewis, 2020). If R j > 1 , then one individual on patch j is producing more than one individual over the entire population and so patch j is defined as a source. To calculate persistence we can use the basic reproduction number, R 0 , which can be calculated as where ρ() is the spectral radius. If R 0 > 1 , then the metapopulation will persist and if R 0 < 1 then the metapopulation will go extinct, a relationship that holds under any of the model formulations considered here (Cushing & Yicang, 1994;Harrington & Lewis, 2020;Li & Schneider, 2002;van den Driessche & Watmough, 2002). The only conditions required are that K be irreducible, which is biologically satisfied if there is some small positive probability that larvae leaving one patch can eventually arrive on any other patch, and that K be nonnegative, which is always satisfied as the entries are numbers of individuals and cannot be negative.
There are two biologically reasonable properties that also exist mathematically under this framework. First, if the population on any single habitat patch can persist on its own, so that k ii > 1, then the entire metapopulation will persist and R 0 > 1 (Harrington & Lewis, 2020). Second, a metapopulation consisting only of sink patches cannot persist and a metapopulation consisting only of sources cannot go extinct. The mathematical underpinning of these relationships is that the spectral radius must be between the minimum and maximum column sums of a matrix, so in terms of our metapopulation quantities It should be noted that both of these properties exist under a deterministic modeling framework; when stochasticity is added neither the persistence nor extinction results are guaranteed.
Having defined patch contribution and persistence in terms of the next-generation matrix, we now demonstrate how to calculate the next-generation matrix for a model with age-dependent demography, as this is the modeling structure commonly used for sea lice on salmon farms. We also present the construction of the nextgeneration matrix for ordinary differential equation models and discrete time models in Appendices 1 and 2, so that the details are present for the most commonly used population models.

| Calculating the next-generation matrix for models with age-dependent demography
Here, we calculate the next-generation matrix for models which allow for the maturation, survival, birth, and dispersal rates to depend not only on the stage and patch location of an individual, but also on the time that they have spent in a stage. This time is often referred to as stage age and so in these models there are two time variables: the global time of the system, t, and the time that an individual has spent in a particular stage, their stage-age, a.
These models can be specified either as McKendrick von-Foerster partial differential equations, integrodifferential equations, or renewal equations (Feng & Thieme, 2000). In all cases, the dependence of the maturation rate on time spent in a stage allows for the addition of more realistic maturation functions where most individuals mature at some intermediate stage-age, or after some minimum time spent in the stage. In contrast, if models are formulated using ordinary differential equations the time in a stage is always exponentially distributed.
However, the specification of the model equations for these models can be rather complicated and so here we construct the next-generation matrix directly from the maturation, survival, birth, and dispersal functions, but the full model and the derivation of the next-generation matrix can be found in Appendix 3 as well as in Harrington and Lewis (2020). Essentially the construction involves tracking the probability that an individual survives through the different stages on a specific patch, the number of larvae that they produce, and the probability that the larvae successfully disperse from one patch from another.
If the maturation rate from stage n to n + 1 is m j n (a), then the probability that an individual has not yet matured from stage n to n + 1 at stage-age a can by given by M j n (a) = exp − ∫ a 0 m j n ( )d . Then, if the probability that an individual survives up to stage-age a in stage n is given by S j n (a), the total probability that an individual is still alive and has not yet matured at age a is S j n (a)M j n (a). To calculate the probability that an individual starting in stage n eventually reaches stage n + 1 , we multiply the rate at which they are maturing at stage-age a, m j n (a), by the probability that they are still alive and have not yet matured at a, S j n (a)M j n (a), and integrate over all stage-ages at which they could have left the stage: If the rate of maturation is constant so that m j n (a) = m j n and the mortality rate is constant, so that S j n (a) = e − j n a , this reduces to which is the case if the model is formulated as a system of ordinary differential equations.
In the terminal stage, if individuals at stage-age a produce larvae at rate b j (a) and the probability that they survive up to a is S j f (a), then to calculate the overall number of larvae that an individual produces over their lifetime we multiply b j (a) by S f (a) and integrate over all a: Therefore, the entries of the next-generation matrix, k ij , which are the total number of new individuals produced on patch i over the lifetime of an initial individual on patch j, will be given by the probability that an initial individual survives to the terminal stage on patch j, multiplied by the total number of larvae it produces, multiplied by the probability that the larvae successfully disperse from patch j to patch i:

| Application: Sea lice on salmon farms in the Broughton archipelago
In this section, we construct a next-generation matrix to determine the contribution of a single salmon farm to the spread of sea lice in a salmon farming region on the west coast of British Columbia, the (5) F I G U R E 2 Map of the 20 historically active salmon farms in the Broughton Archipelago, for which the nextgeneration matrix is calculated.
Broughton Archipelago (Figure 2). To determine the level of sea lice dispersal away from salmon farms, we use a hydrodynamic model, coupled with a particle tracking model, to track sea lice particles released from 20 historical farms in the Broughton Archipelago (Cantrell et al., 2018). The coupled hydrodynamic and particle tracking models have been validated with sea louse counts on salmon farms in the Broughton (Cantrell et al., 2021). We construct one next-generation matrix for the 20 historical farms in the region for which the hydrodynamic particle tracking model was run (Figure 3a), as well as one for the 11 remaining farms in the area after 2023 ( Figure 3b), subject to First Nations and governmental approval (Brownsey & Chamberlain, 2018

| Dispersal
The probability that a sea louse larvae leaves from one farm and successfully arrives on another depends on several factors including ocean current, temperature, and salinity. To accurately capture this probability, it is necessary to use a computational hydrodynamic model that can track the spread of larvae originating from a given farm as well as the dependence of larval survival on temperature and salinity. In order to determine the probability of larvae dispersing between farms we use connectivity matrices from Cantrell et al. (2018). Details on these connectivity matrices can be found in Cantrell et al. (2018). Briefly, these connectivity matrices are calculated by applying Kernel Density Estimation (KDE) to particle tracking simulations to calculate the infectious density of sea lice at each farm, originating from a given farm. The particle tracking simulations are run on the output generated by a Finite Volume Community Ocean Model (FVCOM) which uses data on tides, wind surface heating, and river discharge to simulate three-dimensional ocean velocity, temperature, and salinity (Foreman et al., 2009). In the particle tracking simulation, the survival of sea louse particles is dependent on temperature and their maturation from non-infectious to infectious larvae is dependent on temperature.
The infectious densities of sea lice are calculated for each particle release day by taking daily snapshots of the particle locations of infectious lice for 11 days post-release and then applying Kernel Density Estimation to the accumulated daily snapshots. A connectivity matrix is then calculated for each particle release day, where the entry in row i and column j of the connectivity matrix is the infectious density of larvae over farm i, produced by larvae initially leaving farm j. In Cantrell et al. (2018) the infectious densities were calculated from 24 h of particle releases, where 50 particles were released each hour and so to calculate the infectious density of one initial release particle we divide the entries in each connectivity matrix by 1200 (50 × 24). Then, to create a single connectivity matrix, C, for the 20 farms in the Broughton Archipelago, we take the average over all the connectivity matrices created for particles released between March 14 and July 20, 2009.
The necessary quantity to construct the next-generation matrix is p ij , the probability that larvae leaving farm j will successfully attach on farm i. To estimate p ij from the entries of the connectivity matrix, c ij , there are several assumptions that need to be made. If we assume F I G U R E 3 (a) The next-generation matrix for the 20 historically active farms in the Broughton Archipelago, and (b) the next-generation matrix containing only the farms remaining in the Broughton Archipelago after 2023, subject to First Nations and governmental approval (Brownsey & Chamberlain, 2018). The entries of the next-generations matrices, k ij are the number of new chalimus stage lice produced on farm i from one initial chalimus on farm j. The column sums, R j , are the total number of chalimus produced on all farms from an initial chalimus on farm j and are shown below each column. Likewise, the row sums are the number of new chalimus received by each farm from all other farms and are shown on the left of each row. These numbers should be taken as relative, rather than absolute, as we do not have a very accurate estimate for the arrival rate of sea lice onto farms, β.   1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18  where p j h (x, t) is the two-dimensional density of infectious lice produced from farm j that are still alive in the water column at position x and time t, and β is the arrival rate of lice moving over the farm arriving onto the farm. However, the entries of the connectivity matrix, c ij , are the infectious density of larvae over farm i, produced by larvae leaving farm j. The infectious densities are calculated by applying Kernel Density Estimation to daily snapshots of infectious particles over 11 days, starting at time t = 0, so roughly We can, therefore, roughly calculate p ij from c ij by assuming that the integrals over time and space can be approximated using their Riemann sums, and that the area of a farm is roughly 0.01 km 2 , where Therefore, to calculate p ij from c ij we need to estimate the arrival rate β. However, very little is known about the arrival rate of lice dispersing from one farm to another, and so the estimate presented here is very uncertain and could be orders of magnitude off from the true arrival rate. We assume here that β = 100/day, and thus the It can be seen from Equations (6) and (9) that each entry in the next-generation matrix, K, is a linear function of the arrival rate β, and therefore changing β will not affect the relative ordering of the contributions of different salmon farms, given by R j (Equation 2), it will only affect the absolute magnitude. The basic reproduction rate, R 0 = ρ(K), is also a linear function of β and therefore any increase or decrease in β will result in the same proportional increase or decrease in R 0 . However, the purpose of this application is not to accurately estimate the basic reproduction number, R 0 , for the Broughton Archipelago, or to accurately estimate the contribution of an individual farm, R j , but rather to compare the relative contributions of different salmon farms in the system, and to investigate the effect of environmental variables on the basic reproduction number. Therefore, we present our estimate of the arrival rate for these purposes only, and our value of R 0 found in the results should not be taken as an accurate estimate.

| Demography
Once infectious sea lice larvae attach to their salmonid hosts they must survive and mature through several attached life stages before they can produce offspring. These demographic rates are dependent on salinity and temperature, and thus some salmon farms may be more productive than others due to favorable environmental conditions. To capture the dependence of demography on salinity and temperature we simplify the attached sea lice life cycle down to three main stages: chalimus, pre-adult, and adult. Survival in each stage is salinity-dependent, maturation is temperature-dependent, and egg viability and production depend on both salinity and temperature. The demographic functions that we use are from models which have previously been fit to sea louse population data and are shown in Table 1.
We calculate the on-patch component of the elements of nextgeneration matrix, k ij , by integrating the demographic functions over all time, and refer to this as the productivty of patch j: so that k ij = (Productivity of patch j) × p ij , where p ij is the probability that larvae successfully disperse from patch j to patch i. To determine the specific temperature and salinity-dependent demographic rates at each farm we find the temperature and salinity that each sea louse particle experiences in the particle tracking simulation when initially released from a farm. We then use the average temperature and salinity of particles over all releases.
We also investigate the effect of varying temperature and salinity on the relative growth and persistence of the metapopulation. To do so, we must calculate survival and maturation rates for temperatures and salinities that farms may not experience in the period for which the FVCOM was run. To keep the variability of temperatures and salinities that exists between farms, we multiply the new temperature or salinity at which we want to evaluate persistence, by the ratio of the mean farm temperature or salinity divided by the mean total temperature or salinity experienced by all farms.

| RE SULTS
In this section, we present the next-generation matrix for sea lice populations on salmon farms in the Broughton Archipelago, the construction of which is detailed in Section 2.4, and determine the relative patch contribution of each farm. We use this system to highlight to potential differences between the connectivity matrix, which only contains information surrounding the probability of dispersal from one farm to another, and the next-generation matrix, which combines dispersal between farms and local productivity of sea lice on a salmon farm. We then demonstrate how the next-generation matrix can be used to investigate the effect of changing demographic rates on growth and persistence in this system. Finally, in the context of salmon farm removal from the Broughton Archipelago, we investigate how the removal of habitat patches affects patch contribution and persistence in this sea louse metapopulation. To better understand the details and construction of the nextgeneration matrix, we also present the connectivity matrix for this system in Figure 4 and the productivity (total number of new larvae produced from one attached chalimus louse) of each farm in Table 2.
The (i, j)th entry of the connectivity matrix, c ij , is the infectious density over farm i of lice leaving farm j. The (i, j)th entry of the nextgeneration matrix, k ij (Equation 6), is constructed by multiplying the productivity of farm j (Equation 10) by p ij = β × 0.01 × c ij , the probability that a larvae leaving farm j attaches on farm i (Equation 9). The farms with the largest column sums of the connectivity matrix are, in declining order, farms 18, 2, and 17. However, farm 2 has a higher productivity than farms 17 or 18, and when the productivity of farm 2 is multiplied by the connectivity then it becomes the largest source of sea lice in the region, as identified by the next-generation matrix.
A further look into the connectivity matrix and productivity table provides more insight into the underlying drivers of the contribution of each farm to the sea louse population. Many of the farms with low connectivity also have low productivity, which may be due to the fact that temperature and salinity affect on farm demographic rates as well as survival and maturation in the particle tracking simulation which underlies the connectivity matrix. However, there are certain farms, such as 11 and 15, which have a comparably high productivity compared to their connectivity. These farms are located in favorable environments with respect to temperature and salinity, but low connectivity due to either distance from other farms or unfavorable currents prevents these farms from acting as larger sources of sea lice.
We also examine how temperature and salinity affect the overall growth and relative persistence of the sea louse metapopulation, as shown in Figure 5. The persistence of the metapopulation is determined by the basic reproduction number, R 0 , which is calculated as the spectral radius of the next-generation matrix. As we do not have an accurate estimate for β, we examine the effect of temperature and salinity on the relative change in persistence or growth of the metapopulation, but refrain from commenting on the absolute growth, as measured by R 0 . We can see that as both temperature and salinity increases, the overall growth of the metapopulation increases, and that salinity has a larger effect on metapopulation growth than temperature. What is also interesting, but cannot be seen from the figure, is that as salinity increases, the farm that receives the most lice switches from farm 3 to farm 7.

In light of the removal of salmon farms in the Broughton
Archipelago we also create a next-generation matrix consisting only of the farms which will remain after 2023 subject to First Nations and government approval, shown in Figure 3b, and examine the differences between this matrix and the next-generation matrix with all farms. Most of the farms which are acting as the largest sources TA B L E 1 The maturation, survival, and birth functions used to create the next-generation matrix for sea louse populations on salmon farms in the Broughton Archipelago. The sea louse life cycle is simplified to three attached stages: chalimus, pre-adult, and adult. For the maturation rate out of stage i, i ∈ (c, p) , where c refers to the chalimus stage and p refers to the pre-adult stage.

Description Function Units References
Survival probability of all stages (original next-generation matrix) to R 0 = 2.25. Again these numbers are calculated using a very rough estimate of the arrival rate onto farms, detailed in Section 2.4.2, and thus it is their relative similarity that is important, rather than the absolute magnitude.

| DISCUSS ION
In this paper, we demonstrated how to use the next-generation matrix to calculate the contribution of each habitat patch to the metapopulation and measure the overall persistence of a metapopulation. We detailed the construction of the next-generation matrix under different model structures to demonstrate the breadth of the approach to several systems. We then constructed the next-generation matrix under an age-dependent modeling framework for sea lice populations on

TA B L E 2
The number of new larvae produced on each farm by a single louse starting in the chalimus stage. The first and third rows contain the farm numbers and the second row and fourth rows contain the number of larvae produced.

F I G U R E 5
The effect of temperature and salinity on the overall growth or persistence of the original sea lice metapopulation of 20 farms, as described by the basic reproduction number, R 0 . We do not have a good estimate for the arrival rate of sea lice onto farms, β, and so the R 0 values should only be interpreted relative to each other, rather than as absolute values. salmon farms in the Broughton Archipelago to illustrate how this approach can be applied to a real system. We determined which salmon farms may be acting as the largest sources of sea lice in this region, how the metapopulation will change once certain farms are removed, and examined the effect of temperature and salinity on the relative growth and persistence of this metapopulation.
Next-generation matrices have been used extensively in epidemiology to study the spread of infectious diseases but have recently been introduced in spatial ecology (Huang & Lewis, 2015; Krkošek & Lewis, 2010;Mckenzie et al., 2012) and evolutionary analysis (Hurford et al., 2010). One of the key benefits of using nextgeneration matrices in epidemiology is that the basic reproduction number, R 0 , for a disease can be calculated as the spectral radius of the next-generation matrix, which is often more analytically tractable than calculating the dominant eigenvalue of the full system to determine the spread. In ecology, one main advantage of this approach is that the mathematical calculation of R 0 can be broken down into biologically relevant quantities, for example, the contribution of different dispersal pathways to growth in a population (De-Camino-Beck & Lewis, 2007) or the contribution of populations on different habitat patches (Harrington & Lewis, 2020). While not novel, we hope next-generation matrices can be used more frequently as a simple and easily biologically interpretable method to measure the contribution of local habitat patches to a metapopulation and determine overall persistence.
We are also by no means the first to attempt to calculate the contribution of a local population, classify patches into sources and sinks, or attempt to measure the persistence of metapopulations.
In the context of low densities, Pulliam (1988) defined a source as a habitat patch that would grow in the absence of immigration and emigration and a sink as a habitat patch that would decline in the absence of immigration and emigration. This is similar to only using the entries along the diagonal of the next-generation matrix to classify sources or sinks, except growth or decline was measured after one time step in a discrete time model, rather than one generation.
However, as discussed at the end of Section 2.2, it is possible to have a metapopulation composed only of sinks based on this definition (k ii < 1 for all i), that persists.
Recognizing that dispersal between patches should also be considered when classifying habitat patches as sources or sinks, both Runge et al. (2006) and Figueira and Crowder (2006) defined new metrics to classify habitat patches that track the contribution of adults on a patch in one time step to the total population on all patches in the next timestep. These metrics are similar to our patch contribution metric R j , except they measure the contribution over one time step rather than one generation, similar to using the dominant eigenvalue of the projection matrix A to determine the stability of the discrete system n(t + 1) = An(t), rather than the spectral radius of the next-generation matrix K. However, the calculations can become complicated if the population is stage-structured and the number of stages and/or patches is large (Appendix A, Runge et al. (2006)).
There are other measures of persistence in metapopulations which also track the number of new individuals contributed to the metapopulation after one generation from an initial individual on one patch. Krkošek and Lewis (2010) define a next-generation operator for general heterogeneous populations which tracks the number of new adults produced in the population from one initial adult after one generation. If b j is the reproductive output on patch j, p ij , is the probability of larvae dispersing between patch j successfully arrives on patch i, and a i is the survival to maturity on patch i, then the contribution of patch j to patch i according to Krkošek and Lewis (2010) can be calculated as b j p ij a i . This patch contribution metric cannot be calculated directly from the next-generation matrices shown in this paper unless there is only one stage. Burgess et al. (2014), following Hastings and Botsford (2006), track the number of new recruits on all patches produced over the lifetime of an initial recruit to determine metapopulation persistence in discrete time models. If recruits are defined to be individuals that have survived to maturity, then the entries of their "realized connectivity matrix" are the same as from Krkošek and Lewis (2010) A be a diagonal matrix with entries a j , where p ij , b j, and a j are the same as the preceding paragraph. If we measure generational output starting at the first attached stage, as we do in this paper, the nextgeneration matrix can be written as PAB, if we measure generational output starting at the larval stage, the matrix can be written as ABP, and if we measure generational output starting at the adult stage according to Krkošek and Lewis (2010) or Burgess et al. (2014) then the matrix can be written as APB. This is because when P is multiplied by a diagonal matrix on the right, the entries of the diagonal matrix multiply each column of P and when P is multiplied on the left, the entries multiply each row. Now for any two matrices X and Y, XY and YX have the same eigenvalues, and because matrix multiplication is associative each of the matrices PAB, ABP, and APB all have the same eigenvalues as well, and, therefore, also the same spectral radius.
While the spectral radius of the next-generation matrix determines the overall persistence of the metapopulation and the column sums determine whether a habitat patch is a source or a sink, it is important to note that these quantities (R j and R 0 ) are measured in generational time, rather than chronological time. Therefore, if there are large differences in generation time between patches, there could be situations in which one patch is producing more new individuals per generation than another patch, even if it is producing fewer individuals in chronological time. This may occur if one patch has a higher productivity, but longer generation time, than another patch.
One method to estimate whether patch contribution is the same on generation and chronological time is to estimate the growth rate of the patch, r j , by assuming that at low population densities, the population will grow exponentially and that R j individuals will be produced after one generation, so e r j T j = R j , where T j is the generation time of the patch. If the relative ordering of the growth rates, r j , is the same as the relative ordering of the column sums of the next-generation matrix, R j , then the relative contributions of the habitat patches will be the same in generational and chronological time. If the relative ordering is different, then the classification of habitat patches into sources and sinks using the column sums of the next-generation matrix still applies, only in this case one source patch could be producing more larvae per generation than another source patch, but less in chronological time.
In the case where there are large differences in the number of larvae produced from different patches when measured in generational and chronological time, and the ordering of source or sink patches is important, it may be more suitable to simulate the metapopulation model directly to understand which patches are contributing the most to the growth or decline of the metapopulation. If the metapopulation model is formulated in discrete time, or in continuous time with no age structure, then this is straightforward assuming the number of stages and patches is not too large. In discrete time, patch contribution could also be calculated directly following Runge et al. (2006) and Figueira and Crowder (2006). However, if it is necessary to capture age-dependent maturation and survival, as in the case of sea lice, then simulating the metapopulation model requires either simulating many coupled partial differential equations or converting the model into discrete time. If converting to discrete time, then the number of equations is given by # stages × # patches × age in the longest stage. In the case of sea lice in the Broughton where some stages can last 100 days, this would mean simulating a system of 100 × 4 × 20 = 8000 equations. Additionally, to simulate the model accurately, the time of dispersal between patches must be estimated.
In the case where either the classification of habitat patches into sources or sinks is of primary importance, or when the relative ordering of patch contribution is the same in generational or chronological time, then using the next-generation matrix to determine sources and sinks has many advantages over simulating the metapopulation model. The entries of the next-generation matrix have a direct, simple, biological interpretation in terms of number of new individuals produced after one generation and the framework is the same for discrete time, continuous time, and age-structured systems of equations. The next-generation matrix is easy to calculate, even when maturation and survival in each stage are dependent on stage age. Perhaps most importantly, there may be cases where quantities required to calculate the next-generation matrix are known but the individual rates required to calculate the fulltime-dependent system are not. For example, if the overall lifetime productivity of a patch is known, but not the stage-specific vital rates, or if the overall probability that larvae successfully disperse between patches has been estimated, but not the timing of dispersal, which is the case for sea lice populations on salmon farms in the Broughton Archipelago.
For sea lice on salmon farms in the Broughton, the relative ordering of the population growth rate of each patch, r j , and the column sums of the next-generation matrix, R j , is the same for all but three pairs of farms: farms 2 and 17, 11 and 16, and 4 and 19. The difference in the calculated growth rate r j was <0.0008 between each pair and most farms treat their salmon yearly to reduce sea lice Peacock et al., 2013), so we expect the differences in population size between pairs to be small. Here we used a genera- da to calculate r j . If we were to split the farms into categories, where large sources are those with R j > 3, small sources are those with 1 < R j < 3, and 'sinks' are those with R j < 1, then no farms would switch categories if the rankings were instead based on growth rate. Therefore while a few slight differences in ordering may exist if the output is measured on a chronological scale, the overall findings of which farms are the largest sources of sea lice in the Broughton remain the same irrespective of the timescale used for measurement.
When calculating the next-generation matrix for sea lice populations on salmon farms in the Broughton Archipelago, we believe there are two main sources of uncertainty: uncertainty in the arrival rate, β, and uncertainty in the number of larval sea lice that remain inside the salmon farm in the particle tracking model. In lab studies, the proportion of infective copepodid lice that attach to salmon varies enormously, from as low as 0.5% to as high as >80%, with recent estimates between 12% and 56% depending on temperature (Skern-Mauritzen et al., 2020). Moreover, the time during which copepodids are given to attach in these experiments also varies, from a few seconds, to 8 h. If we assume an exponential waiting time model, P = 1 -e −βt , where P is the proportion of copepodids that have attached to salmon, then these lab studies give estimates of β between 3/day and 47/day, which are lower than the estimate of 100/day which we use in this paper. Any estimates less than 44/day would result in a basic reproduction number, R 0 < 1, which would mean that the overall sea louse population is not increasing.
However, there is abundant evidence that in the Broughton Archipelago sea lice populations on many salmon farms will grow exponentially in the absence of treatment by farm managers Cantrell et al., 2021;Godwin et al., 2021;Krkošek et al., 2012;Rogers et al., 2013). This discrepancy may be due to the fact that while the hydrodynamic and particle tracking models give good estimates of sea louse dispersal between farms (Cantrell et al., 2021), they likely underestimate the number of sea louse larvae that remain inside a salmon farm. In the hydrodynamic model, the location of the farm does not alter the local hydrodynamics that occurs within the farm, but in reality, the presence of salmon swimming within the net pens may alter the local hydrodynamics sufficiently that some proportion of larvae remains within the farm.
To correct this, it would be possible to add a diagonal matrix to the next-generation matrix presented in this paper, where each entry on the diagonal is the productivity of farm j, multiplied by the number of larvae which remain within and successfully attach to salmon on the farm. Because the most productive farms are also the most highly connected in the Broughton, this would likely not change the relative contribution of each farm, but would increase the absolute contribution of each farm as well as R 0 .
In addition to the sources of uncertainty in the next-generation matrix, there are some other technical aspects which should be considered when applying the next-generation matrix to calculate patch contribution of metapopulation persistence in specific systems. In order to use the next-generation matrix to calculate patch contribution or persistence, we are assuming that our system is autonomous and that the demographic rates do not change with time.
In reality for most systems, including sea lice on salmon farms, environmental variables will fluctuate over time, potentially changing the demographic rates of the population (Cantrell et al., 2020). In our case temperature and salinity change over the course of the spring, but we calculate the next-generation matrix using the mean temperature and salinity that sea lice experience during the particle tracking simulation window. Therefore the entries in our nextgeneration matrix may be slightly different than the true number of newly attached lice produced on other farms from one initially attached louse, depending on the exact time that the louse began its lifecycle during the spring. For temporally oscillating systems it is possible to correct for this difference, though the entries of the next-generation matrix no longer have a simple form and must be computed computationally (Rittenhouse et al., 2016).
Over time periods where temperature and salinity are relatively constant, we can use the next-generation matrix to examine the overall effect of environmental variables on the growth of the metapopulation. The demographic rates at each stage depend explicitly on temperature and salinity but alone, or in a full system of equations, it can be difficult to examine the overall effect of changing environmental variables. However, the basic reproduction number R 0 , calculated from the spectral radius of the next-generation matrix, provides a useful metric of the overall effect. We can infer how the growth of the metapopulation may change among seasons, or as the ocean warms. For sea lice on salmon farms in the Broughton Archipelago, the effect of temperature and salinity on R 0 is very similar to previous results found for a single farm (Rittenhouse et al., 2016). With updated temperature and salinity at each farm, we could calculate the change in growth over years in the springtime, which may help explain the recent sea louse outbreaks during warm years in the Broughton Archipelago . When examining the effect of temperature and salinity on R 0 we do not rerun the hydrodynamic and particle tracking models to recreate the connectivity matrices under new temperature and salinity scenarios, as this is very computationally intensive, but we expect connectivity to increase as temperature and salinity increase due to higher survival of sea lice and faster maturation. However, we believe it would be valuable to rerun hydrodynamic models under different projected ocean scenarios, to investigate the precise changes in connectivity that may occur.
Specific to the management of sea lice populations on salmon farms in the Broughton Archipelago, there are several insights to be gained from our results. The first is that the farms in the most productive environments are also the most highly connected, and thus become the largest contributors of sea lice to this sea louse metapopulation. They occur in two main clusters (shown in Figure 3) and both clusters of farms will remain in the Broughton Archipelago in the current removal plan, subject to First Nations and governmental approval (Brownsey & Chamberlain, 2018). It should be noted that these farms may not necessarily be producing the most number of lice compared to other farms in the region at a given time if their louse population is currently lower than other farms, rather they are the farms that have the largest potential to contribute to spread when sea louse numbers are even across farms. However, due to the highly connected nature of these clusters, coordinated treatment between farms in the clusters or all farms in the region could reduce the number of treatments required and number of sea lice produced on all farms .

ACK N OWLED G M ENTS
The authors would like to thank the members of the Lewis lab for many helpful discussions and suggestions. PDH gratefully acknowledges an NSERC-CGSM scholarship, Queen Elizabeth II scholarship, and an Alberta Graduate Excellence Scholarship; and MAL gratefully acknowledges the Canada Research Chair program and an NSERC Discovery Grant.

CO N FLI C T O F I NTER E S T S TATEM ENT
The authors declare no conflict of interest.

DATA AVA I L A B I L I T Y S TAT E M E N T
The data and code used to conduct the analyses reported in this paper are available through the Dryad Digital Repository: https://doi.org/10.5061/dryad.jm63x sjfv (Harrington et al., 2023).

A PPEN D I X 1
Calculating the next-generation matrix for differential equation models Here, we briefly formalize the calculation of the next-generation matrix for differential equation models. We are considering a model for a species with m sessile stages on l patches with a larval stage that can disperse between patches and all other stages are sedentary and confined a patch. Let n j k be the population in sessile stage k on patch j, then the population dynamics at low population densities can be described by where n is a population vector of size m × l describing the population size on all patches at all stages and A is a matrix describing the population dynamics at low population densities. We model the population size at low densities because we are interested in identifying which patches are acting as sources and supporting the population at low densities and which are acting as sinks. Therefore, Equation (11) should be thought of as the linearization about the zero equilibrium of a potentially more complex model that may include density dependence.
The elements of A contain the maturation and death rates of sedentary stages confined to patches as well as the rates of larval dispersal and attachment between patches.
The next-generation matrix for this system can be calculated by first decomposing A = F − V where F is a non-negative matrix with entries that describe the rates of larval birth and probability of dispersal between patches and attachment as the first sedentary phase, and V is a non-singular M matrix (Berman & Plemmons, 1994) with entries that describe the maturation rates between stages and death rates in a stage van den Driessche and Watmough (2002). Because V is a non-singular M matrix, V −1 is non-negative. Following van den Driessche and Watmough (2002) (11) is arranged so that the populations of stage 1 individuals are in the first l rows, then K L will have a l × l submatrix in the upper left-hand corner (Diekmann et al., 2010). This submatrix is the next-generation matrix K, where the elements of K, k ij , give the number of new (stage 1) individuals produced on patch i from one initial stage 1 individual on patch j. In Equation (11), if the maturation rates from stage k to k + 1 on patch j are m j k , the death rates in stage k are d j k , the birth rate in the last stage on patch j is b j , and the probability of successful attachment on patch i as a larvae leaving patch j is p ij , then the entries of K can be written as:

PATCH METAPOPULATION
Here, we illustrate the construction of the next-generation matrix using a two patch example for a species which has two sessile stages. The population of sessile stage k on path j is given by n j k and the maturation rates, death rates, birth rate, probability of dispersal are defined as in the preceding paragraph. The dynamics of the sys-    .
The contribution of patch 1 is, therefore, and the contribution of patch 2 is In order for a patch to be a source (R j > 1), both the on-patch productivity and the connectivity to other patches need to multiply to be larger than 1. If either is too low, that is, if a patch is highly productive but not well connected, or highly connected but not productive, then the patch will be a sink (R j < 1). There are two ways for the entire metapopulation to persist. Either a single patch can persist on its own (k ii > 1), or there must be sufficient production and connectivity within patches such that R 0 = ρ(K) > 1.

A PPEN D I X 2
Calculating the next-generation matrix for discrete time models The construction of the next-generation matrix for discrete time models is very similar to that of continuous time, though with some minor differences. The population dynamics can now be described by where again n is a population vector of size m × l (number of stages × number of patches) describing the population size on all patches at all sessile stages, but now A is a population projection matrix with transition and survival probabilities as well as births.
To calculate the next-generation matrix for the discrete time system, we now decompose A = F + T, where again F is a matrix that contains all the birth rates and probabilities of successful of larval dispersal between patches, and now T is a matrix that contains all of the survival probabilities and transition probabilities between stages. Under this decomposition, the next-generation matrix with large domain, K L , can by calculated by If Equation (14) is arranged so that the populations of all stage 1 individuals are in the top l rows, then again the next-generation matrix, K, will be the l × l submatrix in the upper left hand corner of K L . In the discrete time framework if the probability of transitioning from stage k to k + 1 is m j k , the probability of remaining in stage k is s j k (diagonal entries of A), the fecundity in patch j is b j , and the probability that a larvae leaving patch j successfully arrives on patch i is p ij , then the entries of K can be written as

A PPEN D I X 3
Details of the model with age-dependent demography in Section 2.3 The full set of age density equations for the density of individuals in sessile stage k on patch i at time t and age a, n i k (t, a), is where S j k (a) is the probability of survival in stage k on patch j at stageage a, m j k (a) is the maturation rate from stage k to k + 1, b j (a) is the birth rate of larvae on patch j, and the probability of a larva leaving patch j and successfully attaching on patch i is p ij . M j k (a) is the probability that an individual has not yet matured from stage k to k + 1 and can be calculated as M j k (a) = exp − ∫ a 0 m j k ( )d . Additional details of the model equations can be found in Harrington and Lewis (2020), though there the first larval stage is modeled explicitly and so the set of equations is slightly different.
To construct the next-generation matrix for this system, we need to track the number of stage 1 individuals produced on each patch from one initial stage 1 individual on a given patch. We briefly present the idea of the construction here, and more details can be found in Harrington and Lewis (2020).